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ABSTRACT 


Based on histological analyses and field studies, this 
research describes the reproductive ecology of a 
population of Ninia atrata snakes inhabiting an oil 
palm plantation. Furthermore, through a multivariate 
approach, we explored the main drivers of 
reproductive output in N. atrata. Results showed that 
prey abundance and food intake were crucial 
variables contributing to reproductive output. Multiple 
linear regression models showed that neonates had 
high sensitivity (R?255.2996) to extreme changes in 
climate, which was strongly related to slug and snail 
abundance variability and  microhabitat quality. 
Reproductive cycles were markedly different 
between the sexes, being continuous in males and 
cyclical in females. Despite this variation, 
reproductive cycles at the population level were 
seasonal semi-synchronous. Constant recruitment of 
neonates all year, multiple clutches, high mating 
frequency, and continuous sperm production 
characterized the reproductive phenology of N. 
atrata. |n addition, a significant number of 
previtellogenic females presented oviductal sperm as 
well as uterine scars, suggesting a high precocity in 
the species. The main drivers of reproductive output 
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also differed between the sexes. In females, clutch 
size and secondary follicle variability were highly 
related to stomach bolus volume, fat body area, and 
body mass. In males, height of piles of palm leaves 
and body mass, rather than intrinsic reproductive 
traits, were the main drivers of sperm production. 
Nevertheless, in both cases, the relationship 
between body mass, prey abundance, and food 
intake suggests that N. atrata follows the income 
breeding strategy to compensate for reproductive 
costs and to maximize fitness. 


Keywords: Continuous male reproduction; Clutch 
mass; Income breeding; Iteroparity; 
Spermatogenesis; Oogenesis; Reproductiveeffort; 
EI Nifio-Southern Oscillation (ENSO) 


INTRODUCTION 


Since the early efforts of Fitch (1970) to elucidate the 
reproductive cycles of tropical reptiles, research on the 
reproductive phenology of tropical snakes has increased in 
response to the historical disparity between temperate and 
tropical zone studies (Mathies, 2011). In particular, 
researchers in the South American tropics of Brazil and 
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Argentina have expanded our knowledge about the phenology 
and reproductive traits of more than 10 elapid species 
(Almeida-Santos et al., 1998, 2006; Avila et al, 2010; Valdujo 
et al., 2002), 10 viper species (Almeida-Santos et al., 1999, 
2006; Almeida-Santos & Salomao, 2002; Hartmann et al., 
2004; Janeiro-Cinquini, 2004; Leão et al., 2014; Marques et 
al., 2013; Monteiro et al., 2006; Nogueira et al., 2003), nine 
boine species (Bretona & Chiaraviglio, 2003; Chiaraviglio, 
2006; Miranda et al., 2017; Pizzatto, 2005; Pizzatto & 
Marques, 2007; Rivas et al., 2007), and 26 colubrid species 
(Alencar et al., 2012; Ávila et al., 2006; Balestrin & Di- 
Bernardo, 2005; Bizerra et al., 2005; Braz et al., 2014; da 
Costa-Prudente et al., 2014; Dos Santos-Acosta et al., 2006; 
Gaiarsa et al, 2013; Goldberg, 2004b, 2006 Gomes & 
Marques, 2012; Gualdrón-Durán et al., 2019; Hartmann et al., 
2002; Leite et al., 2009; López & Giraudo, 2008; Marques, 
1996; Marques et al., 2009; Marques & Puorto, 1998; Pizzatto 
et al., 2008; Pizzatto & Marques, 2002; Scartozzoni et al., 
2009; Silva & Vadez, 1989; Vitt, 1996). 

The above contributions have shown that the reproductive 
biology of tropical snakes diverges from the uniformly 
seasonal and highly synchronous patterns seen in species 
from temperate zones. Tropical snakes exhibit both seasonal 
and aseasonal reproductive cycles. However, aseasonal 
reproduction is not the rule in the tropics, in fact, truly 
aseasonal (continuous) reproduction by tropical females is 
rare (Brown & Shine, 2006b; Saint Girons & Pfeffer, 1972; 
Shine, 2003). Tropical snake reproduction includes multiple 
phenologies and high variability among populations, as well as 
within individuals of the same population (Mathies, 2011; 
Pizzatto et al., 2008). Most species also have broad seasonal 
reproductive schedules and exhibit intersexual divergence in 
reproductive cycles (e.g., Naja, Bungarus, Calliophis, Pizzatto 
et al., 2008; Saint Girons & Pfeffer, 1972). 

Despite the enormous efforts to understand the reproductive 
biology of snakes, fewer studies have been conducted on 
tropical fossorial and semifossorial snakes than on terrestrial 
and arboreal snakes (Braz et al., 2014). This is due to the 
cryptic behavior, secretive microhabitats, and lower encounter 
rates of fossorial and semifossorial snakes, which makes them 
an elusive research object. Therefore, their natural history and 
ecology remain poorly understood, which limits our 
assessment of their reproductive seasonality, energy 
acquisition timescale, reproductive expenditure, and sexual 
maturation age or size. In Colombia, the semifossorial snake 
Atractus marthae (Meneses-Pelayo & Passos, 2019) has been 
the only snake species with a detailed reproductive study. 
Atractus marthae populations inhabit the cloud forests (>2 400 
m a.s.l.) of the northeastern Andes of Colombia, with females 
reported to have an asynchronous reproductive stage, 
aseasonal and discontinuous reproductive cycle, and single 
clutch per year, whereas males present spermatozoa in testes 
and ducts, as well as hypertrophied of the sexual segment of 
male kidney (SSK) throughout the year (Guladrón-Durán et 
al., 2019). However, given the dearth of reproductive 
information among Colombian snakes, comparison between 
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the reproductive ecology of snakes from highlands and 
lowlands is difficult, limiting our understanding of the general 
reproductive patterns of snakes in the tropical Andes. 

Ninia atrata (Hallowell, 1845) is a semifossorial tropical 
snake widely distributed in South America. It ranges from 
Western Panama, Colombia, Ecuador, Venezuela to Trinidad 
and Tobago (Angarita-Sierra, 2009, 2014, 2015; Ingrasci, 
2011; McCranie & Wilson, 1995; Medina-Rangel, 2015; Mesa- 
Joya, 2015; Rivas et al., 2012). Despite its broad distribution 
and high abundance in disturbed or transformed habitats, its 
reproductive biology has been largely ignored (Angarita- 
Sierra, 2015; Lynch, 2015). Currently, details on clutch size 
and birth-size in northern South American populations have 
only been reported as anecdotal observations (Lancini, 1979; 
Natera-Mumaw et al., 2015; Roze, 1966). Silva & Valdez 
(1989) reported a hatchling period from June to September in 
populations located on the northern hills of Caracas, 
Venezuela. 

Herein, we explored the reproductive biology of N. atrata in 
order to improve our understanding of the reproductive 
features of semifossorial tropical snakes. Our study consisted 
of three main goals: (1) To provide a detailed description of 
female and male annual reproductive cycles, minimum size at 
sexual maturity, mating frequency, recruitment, clutch size, 
and egg features. (2) To analyze the relationship between 
reproductive traits of N. atrata, prey abundance, and food 
intake, as well as variability during extreme EI Nifo-Southern 
Oscillation (ENSO) climatic events. (3) To explore the main 
drivers of reproductive output of the population under study. 


MATERIALS AND METHODS 


Study area and climate variability 

We studied N. atrata snakes inhabiting the oil palm plantation 
(Elaeis guineensis Jacq. 1897) of Palmasol S. A., located at 
Vereda La Castañeda, the municipality of San Martin, 
Department of Meta, Colombia (N3° 31', W73° 32"; Figure 1). 
This locality has a monomodal climate (rainy season from 
April to November and dry season from December to March) 
with an annual rainfall of the 3 070 mm and high temperatures 
year-round (226 °C). According to the Colombian Institute of 
Hydrology, Meteorology, and Environmental Studies (IDEAM 
Spanish acronym), the El Niño ENSO recorded between 2016 
and 2017 was the strongest of the last 20 years (IDEAM, 
2016). Therefore, climatic variability was categorized as good 
and bad climatic years, where good years represent the 
sampling period from August 2014 to December 2015 without 
ENSO effects, and bad years represent the sampling period 
from January 2016 to April 2017 under ENSO effects. 

We monitored the temperature and relative humidity of 
microhabitats using four thermo-hygrometers (EBI20-TH1 
Ebro&). These devices were placed over three consecutive 
days per sampling visit, with three at the base of different piles 
of palm leaves and one at the base of an epiphytic cushion. 
These sites were randomly selected. A fifth thermo- 
hygrometer was placed at a height of 1.5 m on a randomly 


wrYas Ww3'32 Wrst" Wa" 
Geographic Coordinate System 
Dalum: WGS 1384 


N37 


NSZ 


N3'33 


LEGEND 
BR sapi ved 
Oter batch 
— Streams. 


Forest 
> 


N3'30* 











wry3s' W73'S2 wrisur arro 
Figure 1 Study area 
Oil palm plantation (Elaeis guineensis Jacq. 1897) of Palmasol S. A. 
Sampling batches 8, 9, 13, and 15 are in red 


selected palm tree trunk. We configured the thermo- 
hygrometers to record temperature and relative humidity 
readings on a per minute basis, resulting in a total of 138 240 
measurements and 2304 h of sampling effort per thermo- 
hygrometer. 


Sampling and data collection 

Data were obtained during four plantation production batches 
(8, 9, 13, 15; Figure 1), in which piles of oil palm leaves are 
stacked following pruning. We sampled from August 2014 to 
June 2017, spending three days per sampling visit (0232). No 
sampling was conducted in September 2014 or May 2015 due 
to logistic constraints. 

Snake searching was conducted in a total of 4 563 palm leaf 
piles (microhabitat 1) and 304 epiphytic cushions attached to 
the base of palm trees (microhabitat 2) from 0730 h to 1730 h, 
totaling 1 189 h of sampling effort per researcher. Palm leaf 
piles and epiphytic cushions were searched by raking to a 
depth of 5-15 cm over a 3-5 min period, with all hand- 
captured prey items (e.g., snails, earthworms, slugs, and 
leeches) and snakes recorded. We used the height (cm) of the 
palm leaf piles and epiphytic cushions as surrogate estimators 
of microhabitat quality. This is because greater habitat height 
likely provided better refuge for the snakes due to lower 
microclimate variability, larger number of prey, and greater 
security against predators (Angarita-Sierra, 2015; Angarita- 
Sierra & Lozano-Daza, 2019; Lynch, 2015; Weatherhead & 
Madsen, 2006). 

All caught snakes were measured (snout-vent length (SVL); 
tail length (TL) using a measuring tape (+0.1 cm) and 
weighed (Mass) using a Pesola® dynamometer of 50 g (+0.1 


9). We calculated the sexual size dimorphism index (SSD) 
following Gibbons & Lovich (1990) and determined the ratio 
between the SVL of the larger and smaller sex. The ratio is 
defined as positive when females are larger and negative 
when males are larger (with SSD=1.0 when SVL is equal 
between the two sexes) We recorded snake health 
(presence/ absence of tick, mycoses, or injuries) by external 
examination, as well as sex, umbilicus scars, and secondary 
sexual traits in males (development degree of chin tubercles). 
We also recorded the presence of food and reproductive 
condition by palpation and by contrast light on the snake body 
from the dorsal to ventral surface. 

Snakes captured in batch 13 (n=275) were used in the 
mark-recapture experiments to test whether females could 
produce more than one clutch during the same reproductive 
season. Thus, these snakes were branded on their ventral 
scales following the procedures described by Dorcas & 
Willson (2009) using an Aaron Medical Change-A-Tip cautery 
unit Bovie& (Winne et al., 2006). We determined the sex of 
the branded snakes by inserting a blunt probe through the 
cloaca orifice following procedures described by Blanchard & 
Finster (1933). Afterward, all branded individuals were 
released at the same place in which they were captured. 

We transported gravid females and eggs found in the field 
to the lab in individual terraria. The females were provided 
with food and water ad libitum. The captive environment was 
maintained to resemble the conditions of the oil palm 
plantation (photoperiod regime: 12 h light/12 h darkness; 
temperature: 26.21£0.61 °C and 65.64+3.52 h). Once females 
laid their eggs, we collected, measured, inspected, and 
monitored their temperature throughout incubation 
(26.21+0.61 °C; 65.64+3.52 h) using an infrared thermometer 
(CGHM-H13). Finally, when hatchlings were born, we 
measured, individually marked, and released them with their 
mothers at the place where she was captured. 


Reproductive data acquisition 
In total, 150 snakes from batches 8, 9, and 15 were 
euthanized with an injection of lidocaine 296 (C14H22N20) in 
the heart, then fixed with 1096 formalin, and preserved in 7096 
ethanol. All specimens were deposited in the reptile collection 
of the Instituto de Ciencias Naturales (ICN) of the Universidad 
Nacional de Colombia. We determined sex of the collected 
snakes by abdominal dissection and direct gonad observation. 

Afterward, the reproductive tracts were fixed in 1096 
buffered formalin and then used to make histological slides 
following Luna (1968) to determine reproductive stage, as well 
as spermatogenesis and oogenesis cycles. We registered 
microscopic reproductive traits following the procedures 
described by Krohmer et al., (2004), Balestrin & Di-Bernardo 
(2005), and Ramos-Pallares et al. (2015). Sperm abundance 
and time of spermatogenesis were based on counts (three 
replicates) in cross-sections of seminiferous tubules chosen at 
random using a light microscope under 40x magnification 
(Fox, 1952). 

Based on digital pictures and ImageJ v1.52 software 
(Schneider et al, 2012), we obtained the following 
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macroscopic reproductive variables from sacrificed male and 
female snakes: size of hypertrophy SSK, testicular volume 
(mm?), width of distal end of deferent duct (mm), oviductal 
width (mm), number and diameter of previtellogenic and 
vitellogenic follicles (mm), and width and length of oviductal 
eggs (mm). 

We considered sexually mature males as the smallest (SVL) 
male having spermatozoa in their testes, and sexually mature 
females as the smallest female with vitellogenic follicles or 
oviductal eggs. We calculated relative fecundity (RF) and 
relative clutch mass (RCM) following Iverson (1987) and 
Seigel & Fitch (1984), respectively. We classified female 
reproductive condition as previtellogenic (only translucent tiny 
follicles), vitellogenic (yellowish yolky follicles), ovigerous (with 
oviductal eggs), or vitellogenic and ovigerous (with vitellogenic 
follicles and oviductal eggs simultaneously). We employed 
histological cuts to evaluate and validate the macroscopic 
state of the follicles, which allowed us to allocate an unbiased 
reproductive state for both males and females. Finally, we 
used the presence or absence of sperm in the oviducts and 
infundibulum as a surrogate estimator of mating season. 


Statistical analysis 

We assessed differences in the sex ratio among mature males 
and females each month using the G-test and Chi-square (X?) 
test. We compared size at sexual maturity between males and 
females using t-test and assessed the assumptions of 
normality and homoscedasticity based on Kolmogorov- 
Smirnov and Levene tests, respectively (Guisande-Gonzales 
et al., 2014). We evaluated monthly intersexual variations 
(synchrony) in female and male reproductive stages, as well 
as time variation (seasonality) using the x? test and G-test. 
Likewise, we employed one-way analysis of variance 
(ANOVA) to compare the oviduct distal width between female 
reproductive stages. 

Based on multiple regression analyses of 150 collected 
snakes (females=79, males=71), we explored the main drivers 
of reproductive output of the N. atrata population under study. 
We first considered the following variables: SVL (Var 1), body 
mass (Var 2), primary follicle number (Var 3), secondary 
follicle number (Var 4), fat body area (Var 5), stomach volume 
(Var 6), height of piles of palm leaves (Var 9), testicular 
volume (Var 10), SSK width (Var 11), and width of distal end 
of deferent duct (Var 12). We then measured stomach bolus 
and fat bodies through digital pictures with ImageJ v1.52 
(Schneider et al., 2012). We calculated fat body area (mm?) as 
the sum of each small fat body area attached to the digestive 
and reproductive tract of each individual. We estimated testes 
and stomach bolus volume employing the ellipsoid formula 


(equation 1): 
3, [WG] 1 
aE Ñ Sb 


where, V-testes volume estimated or stomach bolus, w-width, 
and /-length. 
We next used clutch size and sperm count as dependent 
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variables, with all remaining variables considered 
independent. Both dependent and independent variables were 
square  roottransformed using the Tukey's staircase 
transformation method described by Erickson & Nosanchuk 
(1977). 

We evaluated assumptions of normality, autocorrelation, 
and homoscedasticity using the Kolmogorov-Smirnov test, 
Durbin-Watson test, and Breusch-Pagan test, respectively. 
We tested multicollinearity between the previously named 
variables using the variance inflation factor (VIF) with a 
threshold of 10. To select the "best" regression model based 
on the variables evaluated, we employed the Akaike 
Information Criterion (AIC; Akaike, 1973). 

We analyzed recruitment variability during ENSO events, 
employing neonatal abundance as the dependent variable and 
height of piles of palm leaves (Var 9) and prey abundance 
observed in microhabitats (snails, slugs, earthworms, leeches; 
Angarita-Sierra & Lozano-Daza, 2019) as independent 
variables. 

Finally, to estimate the contribution of all independent 
variables to the regression models, we employed a 
hierarchical partitioning method (Chevan & Sutherland, 1991). 
All statistical analyses were performed using Rwizard 3.0 
(Guisande-Gonzáles et al, 2014) and the following R 
packages: StatR (Guisande-Gonzáles et al., 2014), hier. part 
(Walsh & MacNally, 2015), nortest (Gross, 2015), Imtest 
(Hothorn et al., 2017), and usdm (Naimi, 2015). 


RESULTS 


From all batches sampled, we caught a total of 425 specimens 
of N. atrata (males-209, females-216), with an SSD index of 
0.16. The sex ratio showed no significant differences between 
months (x?=15.89, P=0.145, n=425), climate season (y7=2.30, 
P=0.1286, n=425), or batches sampled (6470.93, P=0.61, 
n=425). |n contrast, significant differences in snake 
abundance were observed between good and bad years 
(x2=176.15, P<0.0001, n=425; Figure 2). Health condition was 
recorded as a potential variable affecting reproductive output. 
However, through the whole sampling period, only 11 
specimens (2.6%) suffered poor condition (ticks n=1; mycoses 
n=8; injury=2). Therefore, these variables were excluded from 
analyses and we considered the N. atrata population under 
study to be healthy. 


Female reproductive cycle and activity 

According to the macroscopic reproductive traits observed in 
79 N. atrata females, the smallest female with vitellogenic 
follicles, indicating sexual maturity, was 270 mm SVL. All 
females larger than 270 mm SVL (48%) were in some stage of 
reproduction (Table 1). Females were asynchronous in 
reproductive stage between months or climate years (X? 
months=27.374, P=0.44; X? climate years=3.05, P=0.38). All 
reproductive stages were observed throughout the year. 
Oviduct width exhibited significant differences among female 
reproductive stages (ANOVA F=13.7, P<0.000 1), with the 
oviduct being less wide in previtellogenic females than the 


I7 Good years ------------------------- | —————-7————--- Bad years -----------------------—------ - 


Rainy season Dry season Rainy season 





i i 
i H 
' H 
30- | i 
H H 
i i 
H H 
ie H i 
o i H 
x H 1 
o 1 l 
$5 20 ! | 
L L H H 
o H i 
o i | 

g i 

: : | 
v : j 
c " 
3 ; l 
G H l 
L d j 
g 10 } i 
E] d i 
e i I 
9 i 1 
2 M , 
« h h 
d 1 
4 
0- i j 


Aug-2014 





Dec-14 Feb-15 Jun-15 Jul-15 Oct-15 Dec-15 Feb-16 Apr-2016 Jun-16 Aug-2016 Oct-16 Dec-16 
Month 


Dry season 


Rainy season Dry season 


* Males 
= Females 


1 
i 
i 
H 
H 
1 
i 
i 
H 
i 
1 
i 
: 
H 
H 
H 
$ 
H 
H 
i 
i 
H 
i 
H 
i 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
H 
1 


Feb-17 Apr-2017 


Figure 2 Abundance and population sex ratio of Ninia atrata through sampling period (2014-2017) 


Table 1 Mean values of female reproductive features 





emae Oviduct length No. primary No. seconda 
reproductive SVL (mm) Mass (g) g Sperm in oviduct Uterine scar : p y ; y No. eggs 
stage (mm) follicles follicles 
205.05 4.94 4.00 Present (40%) Present (15%) 5.71 
Previtellogenic (120-246) (1.08-9.42) (0.91-9.75) Absent (60%) Absent (85%) (1-10) Absent Absent 
n=38 n=34 n=38 n=43 n=38 n=38 
313.57 12.01 9.19 Present (87.5%) Present (20%) 7.21 2.36 
Vitellogenic (270-353) (8.33-17.09) (2.59-16.85) Absent (12.5%) Absent (80%) (3-10) (1-5) Absent 
n=14 n=14 n=14 n=16 n=10 n=14 n=10 
335.08 15.88 9.06 Present (83.3%) Present (100%) 8.66 3.25 
Ovigerous (278-370) | (13.33-19.8) (4.10-13.49) Absent (16.7%) n=12 K (5-13) Absent (2-4) 
n=12 n=11 n=12 n=12 n=12 n=12 
] : 315.68 15.02 10.13 $ o 10 1.8 2.8 
kier MN (320-370) — (12.35-19.1) (5.19-14.4) oe m i dam (100%) quod (1-3) (2-4) 
9 n=5 n=5 n=5 E ~ n=5 n=5 n=5 
remaining reproductive stages. However, differences in April to November (rainy season), with abundance and size 


oviduct width among females at ovigerous, vitellogenic, and 
vitellogenic-ovigerous stages were not significant (ANOVA F= 
0.14, P=0.86). Thus, an oviduct width larger than 6.14+3.85 
mm indicated sexual maturity. Nonetheless, a high degree of 
oviduct width overlap was observed among female 
reproductive stages, suggesting that this macroscopic 
character is not an accurate predictor of sexual maturity 
(Figure 3). 

In the 79 adult females examined, primary follicles were 
present throughout the year. However, in mid-June, most 
follicles began to enlarge, reaching a maximum in November 
to January until the beginning of the dry season (Figure 4A). 
Likewise, the greatest abundance of primary follicles was 
observed in January to March and September to November. 
Consequently, secondary follicles were only observed from 


increasing gradually, reaching the greatest size in the October 
to November period (Figure 4B). 

We observed mating signals at all female reproductive 
stages. However, the number of immature females exhibiting 
mating signals was significantly lower than that of mature 
females (x?^-20.14, P«0.001, n-79) As expected, the 
frequency of uterine scars was significantly higher in females 
at ovigerous and vitellogenic stages (X?-44.32, P«0.001, 
n=79). Notwithstanding, uterine scars also were observed in 
six previtellogenic females, four of which had sperm inside the 
infundibulum and two of which had sperm inside the oviduct. 
Mature females had sperm inside their infundibulum or oviduct 
almost the entire year (except December), indicating that 
copulation is continuous even for females not ready to mate. 
We observed eggs throughout the year, though the greatest 
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Figure 3 Variation in oviduct width 
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differences between female reproductive stages 


abundance was recorded in September to December 
(Figure 4C). Based on the incubation of three clutches 
oviposited by three gravid females, the estimated birth-time 
was 10841.41 d/egg. Clutch size ranged from 1-4 eggs 
(2.44x1.02, n=34). We recaptured two female snakes who 
produced eggs twice in the same reproductive season (first 
capture in July and recapture in November). Also, during the 
dissection of snakes collected in August to November from 
batches 8 and 9, we recorded five females with vitellogenic 
follicles and oviductal eggs simultaneously, indicating that 
females could produce more than one clutch per reproductive 
season. 

Neonates were also observed throughout most of the year 
(except July), with three abundance peaks during the 
recruitment season. The greatest recruitment peak occurred 
during the early dry season (January-February). A second 
moderate recruitment peak was observed in September, and a 
third pronounced recruitment peak occurred at the end of the 
rainy season (November). The greatest dearth occurred 
during the first half of the rainy season (June—July; Figure 4D). 
Based on the birth of eight neonates, birth-size was estimated 
to be 114.63410.69 mm SVL and 1.91+0.74 g of body mass. 
Likewise, RCM and RF were highly variable and ranged from 
4.32% to 8.54% (7.2141.14%, n=16) and 2.72% to 12.74% 
(7.2843.01%, n=34), respectively. 

We observed a remarkable decrease in neonates between 
good and bad climate years. Based on the presence or 
absence of the ENSO, the number of neonates declined 
significantly, from 57 in good years to four in bad years 
(10.557, P=0.001, n=61). Despite this, we observed a clear 
synchronization between recruitment peaks and prey 
abundance in years without ENSO effects. Increased snail 
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and slug abundance coincided with increased neonate 
abundance over the same time period (Figure 5). This 
relationship was confirmed by multiple regression analysis, 
with neonate abundance being strongly correlated with slug 
and snail abundance, but not with other prey (R?-0.46, P- 
0.032; Figure 6). Nonetheless, when the height of piles of 
palm leaves was included in the analysis, it provided a better 
model, explaining 55.29% (P=0.011) of the neonate 
abundance variability observed (Table 2). 


Male reproductive activity 

Ninia atrata males exhibited early sexual activity. The smallest 
male with sperm in their testes, indicating sexual maturity, was 
145 mm SVL, and the largest male without sperm was 212 
mm SVL. These males represent the extreme body-size limits 
of sexual maturity. However, 98.60% (n=71) of males larger 
than 187 mm SVL had sperm in their testes and deferent duct. 
The smallest male with metamorphosing spermatocytes was 
137 mm SVL. 

Based on macroscopic examination of the male gonads, 
testes showed noticeable size variation throughout the year. 
Testicular size gradually changed, decreasing from February 
to August and increasing from August to November. Maximum 
volume was attained in April (beginning of rainy season) and 
October-November (end of rainy season) and was 136.4% 
greater than the minimum testicular volume observed in 
August (mid rainy season; Figure 7A). Despite a lack of 
mature male samples in December, January, and March, our 
data suggests that testicular volume declines at the beginning 
of the dry season but increases mid-way through. 

Likewise, macroscopic sexual features such as SSK and 
distal end of deferent duct exhibited a similar monthly 
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Figure 4 Monthly variability of female reproductive traits and neonates 
Black dots represent means, black bars represent standard deviation. A: Absolute values of primary follicles. B: Absolute values of secondary 
follicles. C: Absolute values of eggs found during sampling period and in dissected reproductive tracts. D: Absolute values of neonates during 


sampling period 


variability pattern as observed for testicular volume 
(Figure 7C-D). Indeed, monthly variability of these traits was 
closely related to testis size (Rss? =0.68, P«0.000 1, n=71; 
Rueterent duc? = 0.45, «0.000 1, n=71). In contrast, sperm 
production was not correlated with monthly variability in 
macroscopic male sexual features (Rssk? =0.031, P=0.15, 
n=69; Rrestis volume =0.028, P=0.18, n-71; Rueterent duct 70.033, 
P=0.14, n=71), indicating that testicular volume, SSK 
hypertrophy, and deferent duct width are not concordant with 
spermatogenic activity. Sperm production was present for 


most of the year, including the dry and rainy seasons. 
Specifically, production increased gradually from April to 
November and reached a maximum in July-August (mid rainy 
season) without significant decline once production began 
(Figure 7B). 

Conversely, significant differences in the size of males with 
either weak or prominent chin tubercles were found (ANOVA 
F=37.28, P<0.000 1, n=71). The weak condition was generally 
associated with males within the SVL range of 135-241 mm 
(168.16+2.02, n=12), although two larger males (SVL=276 
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and 280) exhibited this condition even though they had 
reached the minimum size of sexual maturity. In contrast, the 
prominent condition was generally associated with males 
within the SVL range of 183—354 mm, in accordance with the 
results that 98.60% (n=71) of males of this size have sperm in 
their testes and deferent duct. Despite this, three smaller- 
sized males (SVL=146, 154, and 175) also had prominent chin 
tubercles. In fact, the SVL ranges of chin tubercle condition 
showed a high degree of overlap (46.53%), indicating that this 
secondary sexual character may not be an accurate predictor 
of male reproductive stage. 


Reproductive cycles and mating 
We found notable differences in the reproductive cycles 
between the sexes. First, males showed a continuous cyclical 
pattern, in which spermatogenesis, gonads, and SSK were 
active throughout the year. Although they did show reduced 
activity during the dry and mid rainy seasons, they never 
displayed total regression or quiescence. In contrast, females 
showed a cyclical pattern in which oogenesis, gonads, and 
accessory organs become inactive or absent during the dry 
season. Second, size at sexual maturity was significantly 
different between the sexes (9.54, P«0.000 1, n=150), with 
males and females attaining sexual maturity at 5696 and 8696 
of mean adult SVL, respectively. Finally, higher sperm and 
vitellogenic follicle production were not synchronized. While 
maximum sperm abundance occurred from July to August 
(mid rainy season), maximum vitellogenic follicle abundance 
occurred from October to November (end of rainy season). 
Despite the divergence in reproductive cycles between 
sexes, sperm production and follicle maturation patterns 
indicated that the reproductive cycle was seasonal at the 
population level. Both sexes presented a synchronized 
increase in reproductive output through the rainy season, with 
highest abundance from June to November. Even though no 
mating behaviors were observed among individuals of N. 
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Figure 6 Multiple regression models related to neonate and prey 
abundances 

Plots depicting linear regressions between neonate abundance (A) 
and prey abundance (B) Relative contributions of variables. Ln 
Neonates: Logarithm-transformed 
Earthworms: Logarithm-transformed earthworm abundance Ln Snails: 
Logarithm-transformed snail abundance. Ln Slugs: Logarithm- 
transformed slug abundance 


neonate abundance. Ln 


atrata, the high frequencies of oviductal sperm, as well as 
sperm production and follicle maturation, suggest two main 
mating pulses, one at the beginning of the rainy season (April) 
and one at the end of the rainy season (October-November). 
However, mating signals were present all year, including the 
dry season (Figure 8). 

The seasonal reproductive cycle of N. atrata was closely 
correlated with the strong climate seasonality of the study 
area, as well as prey and hatchling abundance. For instance, 
the greatest recruitment peak occurred in the mid dry season, 
which coincided with the increase in snail abundance and 
greatest production of primary follicles, whereas the greatest 
recruitment dearth was observed in the mid rainy season, 
which coincided with the decline in snail abundance and 
maximum sperm production (Figures 7—8). 


Main drivers of reproductive output in Ninia atrata 

The main drivers of reproductive output in N. atrata diverged 
between the sexes. In females, multiple regression analysis 
indicated that clutch size was strongly correlated with almost 
all reproductive traits evaluated (Table 2). However, among 


Table 2 Multiple regression analysis models 


Nor. test P Hom.testP Aut. test P 
































Model R? AIC dAIC df 
value value value 

Clutch size versus female reproductive drivers 
Clutch size vs. Var2+Var3+Var4 —90.62 0.0 1 
Clutch size vs. Var2+Var3+Var4+Var5+Var6 —89.16 —1.46 1 
Clutch size vs. Var1+Var2+Var3+Var4+Var5+Var6 0-70 -87.28 -2.98 1 0.07 0.35 0.63 
Clutch size vs. 
Vari IVar Var3* Var4* Varb* VarG* Varg Due. com 1 
Secondary follicles versus female reproductive drivers 
Secondary follicles vs. Var2+Var5+Var6 11.18 0.0 1 
Secondary follicles vs. Var2+Var5+Var6+Var9 11.69 0.51 1 

. 0.42 0.17 0.36 0.87 
Secondary follicles vs. Var2+Var3+Var5+Var6+Var9 13.24 2.06 1 
Secondary follicles vs. Var1+Var2+Var3+Var5+Var6 15.21 4.03 1 
Sperm count versus male reproductive traits and environmental variables 
Sperm count vs. Var2+Var6+Var9 —166.21 0.0 1 
Sperm count vs. Var1+Var2+Var6+Var9 —165.20  -1.01 1 
Sperm count vs. Var1+Var2+Var6+Var9+Var10 Tm —163.32 -2.89 1 i 0:33 We 
Sperm count vs. Var1+Var2+Var6+Var9+Var10+Var11 —161.38 -4.83 1 
Neonates versus prey abundances 
Neonates vs. Slugs+Snails —9.06 0.0 1 
Neonates vs. Slugs+Snails+Leeches 0.46 -7.60 1.46 1 0.34 0.50 0.46 
Neonates vs. Slugs+Snails+Leeches+Earthworms —6.41 2.65 1 
Neonates versus prey abundances and height of piles of palm leaves 
Neonates vs. Var9*Snails -11.59 0.0 1 
Neonates vs. Var9+Snails+ Slugs —11.49 —0.10 1 
Neonates vs. Var9+Snails+ Slugs+Leeches 0.55 —10.55 —1.04 1 0.83 0.43 0.92 
Neonates vs. Var9+Snails+ -8.65 -2.94 


Slugs+Leeches+Earthworms 

AIC: Akaike information criterion, employed to select“best model", was used to test whether environmental factors rather than intrinsic reproductive 
traits are main drivers of reproductive output. Var1: Snout-vent length, Var2: Body mass, Var3: Primary follicles number, Var4: Secondary follicles 
number, Var5: Fat body area, and Var6: Stomach bolus volume, Var9: Height of piles of palm leaves, Var10: Testicular volume, Var11: Width of 


sexual segment of kidney, and Var12: Distal width of deferent duct. 


R?: Proportion of variance for reproductive output explained by 


microenvironment or reproductive intrinsic trait variables. Nor. test: Kolmogorov-Smirnov test for normality; Hom. test: Breusch -Pagan test for 


homoscedasticity; Aut. test: Durbin-Watson test for autocorrelation. 


these variables, the "best model" was comprised of the 
number of primary and secondary follicles and body mass, 
which explained 70.63% (P<0.000 1) of clutch size variability. 
Given the importance of secondary follicles in clutch size, a 
second multiple regression analysis was carried out exploring 
the relationships among secondary follicles, maternal traits, 
and height of piles of palm leaves (Figure 9). As a result, 
secondary follicle variability was highly correlated with 
stomach bolus volume, fat body area, and body mass. These 
variables composed the "best model" and explained 44% (P= 
0.003) of secondary follicle variability. Similarly, female SVL, 
but not the remaining variables, was significantly related to 
egg mass (F=7.64, P=0.014, n=17), indicating that larger 
females produced heavier eggs. 

In males, body mass, height of piles of palm leaves, and 
stomach volume, rather than intrinsic reproductive traits, 
showed the greatest contribution to sperm production (R?- 
0.198, P<0.001; Table 2). This result agrees with the 
discordance observed between sperm production and monthly 


variation in testicular volume and size of SSK (Figure 10). 


DISCUSSION 


In general, the reproductive ecology of N. atrata followed 
typical patterns reported in tropical snakes. First, early male 
maturation at a smaller size than females agrees with the 
common pattern among oviparous, small or median sized 
dipsadid species (Dos Santos-Acosta et al., 2006; Parker & 
Plummer, 1987; Pizzatto et al., 2008). Indeed, Goldberg 
(2004a) reported a similar maturation size for N. maculata 
(Peters, 1861) in which the smallest spermiogenic male was 
179 mm SVL and the smallest vitellogenic female was 190 
mm SVL. 

Second, the reproductive cycles of both males and females 
were asynchronous, whereas the reproductive cycle at the 
population level was seasonal semi-synchronous, which 
agrees with the patterns observed in several tropical snakes 
with diverse phylogenetic histories, such as Atractus marthae 
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Figure 7 Monthly variability of male reproductive traits 


Black dots represent means, black bars represent standard deviation. A: Absolute values of testis volume. B: Average values of sperm count per 


individual. C: Width of sexual segment of male kidney. D: Width of distal end of deferent duct 
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Figure 8 Oogenesis and spermatogenesis cycles in Ninia atrata 
Females: 79, Males: 71. Boxes: Observed, dotted line: Inferred 
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Figure 9 Multiple regression models related to clutch size and secondary follicles with environmental variables and intrinsic 


reproductive traits 


A, B: Plots depicting linear regressions between clutch size and secondary follicles with variables, respectively. C, D: Relative contributions of 
variables composing"best regression model"for clutch size and secondary follicles, respectively. Var2: Body mass; Var3: Primary follicles number; 
Var4: Secondary follicles number; Var5: Fat body area; Var6: Stomach bolus volume 


(Meneses-Pelayo & Passos, 2019), Sibynomorphus spp., 
(Schlegel, 1837), Atractus reticulatus (Boulenger, 1885), 
Drymobius margaritiferus (Schlegel, 1837), Dipsas albifrons 
(Sauvage, 1884), Mastigodryas melanolomus (Cope, 1868), 
Micrurus lemniscatus (Linnaeus, 1758), and others (Goldberg, 
2006; Marques et al., 2013; Pizzatto et al., 2008). 

Third, the presence of several vitellogenic-ovigerous 
females during the reproductive season, constant recruitment 
of neonates year-round, frequent mating evidence all year 
independent of season, and presence of previtellogenic 
females with sperm in their oviduct or infundibulum agrees 
with the patterns reported for numerous tropical snakes, such 
as Erythrolamprus aesculapii (Linnaeus, 1758), Erithrolamprus 
bizona (Jan, 1863), Mastigodryas bifossatus (Raddi, 1820), 
Tropidonophis mairii (Grey, 1841), and others, in which 
multiple clutches, high mating frequency, and continuous 
sperm production characterize their reproductive phenology 
(Brown & Shine, 2006a; Goldberg, 2004b, 2006 Marques, 
1996). This suggests that immature females are willing or are 
forced by males to mate. 

In contrast, the body size-fecundity relationship observed in 
N. atrata moves away from the expected correlation between 
SVL and reproductive output in tropical snakes (Miranda et al., 


2017; Shine & Madsen, 1997). In N. atrata, SVL was shown to 
be a poor morphological predictor of fecundity. In particular, 
SVL was only significantly correlated with egg mass, 
presumably reflecting physical constraints on clutch volume 
(Shine, 1991). In contrast, body mass was shown to be a 
better morphological predictor for both sexes, as this trait was 
persistently selected in all regression models assessed. Male 
body mass had the strongest contribution (25096) to sperm 
production. Similarly, female body mass in all regression 
models evaluated explained 2096 of the reproductive output 
and occupied the second or third place of importance, after 
traits such as secondary follicle number, stomach bolus 
volume, or fat body area. Nonetheless, sperm production was 
poorly explained by the variables assessed (?-19.8894). 
Thus, future studies should clarify which microenvironment 
variables or intrinsic reproductive traits can determine the 
reproductive output observed in N. atrata males. 

Food intake and prey abundance were crucial variables 
contributing to reproductive output in N. atrata. The high 
association of vitellogenic or ovigerous females with stomach 
content, as well as the great importance of stomach bolus 
volume as an explanatory variable of secondary follicle 
variability and sperm production, rather than fat bodies or 
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Figure 10 Multiple regression models related to sperm 
production with environmental variables and intrinsic 


reproductive traits 

A: Plots depicting linear regression models between sperm count with 
variables. B: Relative contributions of variables composing "best 
regression model"for sperm count. Var2: Body mass; Var6: Stomach 
bolus volume; Var9: Height of piles of palm leaves 


SVL, indicated a strong link between N. atrata reproductive 
output and food intake. In addition, the strong link observed 
between recruitment peaks and prey abundance is in 
agreement with the results of Shine (2003), who argued that 
even in the tropics, seasonal reproductive cycles would be 
favored because they reflect the variability of operative 
environmental factors, as well as the temporal shifts in 
reproductive trade-offs. 

Similarly, the remarkable decrease in neonates and 
significant differences in snake abundance during ENSO 
events emphasized the high sensitivity of N. atrata to extreme 
climate changes. The variabilities of operative environmental 
factors such as height of palm leaf piles and prey availability 
(mainly snail abundance) were demonstrated to be the main 
drivers of N. atrata abundance variability (Angarita-Sierra & 
Lozano-Daza, 2019). 

The above evidence suggests that N. atrata populations 
follow an income breeding strategy (Jónsson, 1997) in order to 
compensate for the demands of reproduction and to maximize 
fitness. This result agrees with the observed feeding patterns, 
i.e., largest females had significantly higher stomach bolus 
volumes than the largest males, intersexual dietary 
divergence, and notable disparity in food intake between the 
sexes (Angarita-Sierra & Lozano-Daza, 2019). Moreover, 
unlike tropical snakes that exhibit a seasonal shift in their 
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reproductive strategies (from capital breeding to income 
breeding or vice versa, e. g., Tropidonophis mairii (Brown & 
Shine, 2006a)), N. atrata maintained the same reproductive 
strategy, despite the extreme climatic variability due to ENSO. 

Income breeders are rare among ectotherms because of the 
high energetic costs associated with the muscle and organ 
maintenance needed to sustain their highly active mode of 
feeding (Bonnet et al., 1998). However, the N. atrata 
population in the oil palm plantation had a huge amount of 
prey available year-round, allowing them to reduce their 
energetic costs of acquisition without exceeding the energetic 
costs that would have to be expended in reproduction or other 
activities (Angarita-Sierra & Lozano-Daza, 2019). Hence, the 
income breeding strategy observed in this population could be 
strongly influenced by habitat. However, more empirical data 
are needed to elucidate whether reproductive strategy 
switches with habitat type or is a conservative life-history trait 
of the species. 

Finally, macroscopic reproductive traits were shown to be 
an inaccurate proxy for reproductive activity in both sexes. In 
particular, testicular volume and SSK were not associated with 
sperm production. Likewise, oviduct distal width and presence 
of chin tubercles exhibited wide variance, which made it 
difficult for accurate maturity size determination. The accuracy 
of macroscopic reproductive traits as a proxy for reproductive 
cycle evaluation has been questioned previously (Braz et al., 
2014; Mathies, 2011). In fact, it has been observed that 
histological analysis invalidates macroscopically determined 
maturity in fish, lizards, and fossorial snakes (Boretto & 
Ibargüengoytía, 2006; Braz et al., 2014; Fernandez et al., 
2017; Vitale et al., 2006). Nonetheless, this topic has been 
poorly explored in snake reproductive studies due to the 
historical and widespread use of macroscopic reproductive 
traits in comparisons between different studies (Pizzatto & 
Marques, 2006, 2007). Therefore, for future comparison 
between N. atrata populations or related taxa, we recommend 
employing histological assessment to avoid spurious results 
that could distort the relationships among reproductive cycles, 
environmental factors, and reproductive trade-offs. 
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